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We construct the effective lowest-band Bose-Hubbard model incorporating interaction-induced on-site cor- 
relations. The model is based on ladder operators for local correlated states, which deviate from the usual 
Wannier creation and annihilation, allowing for a systematic construction of the most appropriate single-band 
low-energy description in form of the extended Bose-Hubbard model. A formulation of this model in terms of 
ladder operators not only naturally contains the previously found effective multi-body interactions, but also con- 
tains multi-body induced single particle tunneling, pair tunneling and nearest-neighbor interaction processes of 
higher orders. An alternative description of the same model can be formulated in terms of occupation-dependent 
Bose-Hubbard parameters. These multi-particle effects can be enhanced using Feshbach resonances, leading 
to corrections which are well within experimental reach and of significance to the phase diagram of ultracold 
bosonic atoms in an optical lattice. We analyze the energy reduction mechanism of interacting atoms on a local 
lattice site and show that this cannot be explained only by a spatial broadening of Wannier orbitals on a single 
particle level. 
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I. INTRODUCTION 

Ultracold atoms in optical lattices are an ideal testing 
ground for models in solid-state physics due to the large de- 
gree of control over external and internal parameters of these 
many -body systems [lj. On the one hand, these systems are 
very promising as analog quantum simulators for gaining fur- 
ther understanding of complicated solid state systems [2]-[5]], 
whereas on the other hand completely new models (e.g. with 
further internal degrees of freedom, different quantum statis- 
tics, etc.) can be realized in a very clean and controlled fash- 
ion. Specifically, a large focus has been on ultracold bosonic 
atoms in optical lattices, which are well described by the Bose- 
Hubbard model [2j. The first milestone was the experimen- 
tal observation of the superfluid-Mott insulator transition QQ. 
With the ever increasing precision in recent experiments |7]- 
9j, as well as development of new probing techniques and re- 
markable technical advances H10l - fT3'j , it has become possible to 
observe effects beyond the standard Hubbard model. Specifi- 
cally, a density dependence of the interaction parameter U has 
been observed by using quantum phase revival spectroscopy 
BUI, which has been predicted and described using effective 
many-body interactions. A recent experiment using multi- 
band spectroscopy to investigate the effect of bosons in a Bose- 
Fermi mixture found a significant reduction of the fermionic 
tunneling energy J II 1411 . 

While the Fock space spanned by the Fock states generated 
by the full multi-band single-particle Wannier orbitals is a per- 
fectly valid basis for the interacting many-body system, where 
by construction the parameters are density-independent, it is 
customary to work in an effective single band basis. However, 
such a description requires a density dependence of the pa- 
rameters, or alternatively the introduction of effective higher 
order terms as will be shown. 

It has been proposed II15H17II that the density dependence 
of the bosonic tunneling parameter induced by the Bose-Fermi 
interaction can explain the shift in the bosonic superfluid-Mott 



insulator transition observed in Bose-Fermi mixtures I10II11I . 
This topic is still under debate with an alternative cause sug- 
gested to be the heating of the system as the lattice is ramped 
up II18H20II . Furthermore, several new phases have recently 
been predicted for the effective single-band density-dependent 
Bose-Hubbard model 121112211 . 

In this work we rigorously derive and define the effec- 
tive lowest-band representation used in these previous works, 
where the localized many-body low-energy states are dressed 
with contributions from higher bands, analogous to the 
dressed state basis in quantum optics. We define new lad- 
der operators connecting only states within this dressed low- 
energy manifold, which exactly fulfill bosonic commutation 
relations. For finite interaction strength \g\ > 0, these do 
not coincide with the usual single-particle Wannier creation 
and annihilation operators, as has commonly been assumed 
10 [8] [10] [JJ1 . We give the exact prescription for transform- 
ing between the multi-orbital Wannier and the dressed single- 
band basis in the low-energy regime. This transformation is 
also vital to translate any operator into the new basis, which 
is usually given in the real space, Bloch or Wannier represen- 
tation, e.g. observables, additional terms in the Hamiltonian 
or perturbations. On a local level, our transformation recovers 
the effective multi-body interactions found in [23] and further- 
more allows for a systematic treatment of all non-local terms. 
Here we identify the counterparts of local multi-body inter- 
actions: multi-particle induced tunneling and correlated-pair 
tunneling terms arising from the usual bosonic interacting lat- 
tice Hamiltonian. 

This article is organized as follows: in Sec. jll} we juxta- 
pose the multi-orbital Wannier and the effective single band 
descriptions and introduce the basis states of the latter. In 
Sec. (Hi) and Sec. (TV) we define the low-energy subspace and 
the new effective bosonic ladder operators, from which the 
transformation properties are derived. Subsequently, they are 
applied in the systematic derivation of additional terms to the 
standard Bose-Hubbard model and are shown to give rise to n- 
particle induced single- and correlated two-particle tunneling 
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in in Sec. jV A\ and in Sec. (VB} respectively, as well as multi- 
body nearest neighbor interactions in Sec. (VCj . Finally, we 
investigate the main energy reduction mechanism in Sec. jVII} , 
showing that mutual particle avoidance visible in the second 
order correlation function is more important than the com- 
monly used explanation of broadened single-particle orbitals. 



the nearest neighbor tunneling energy along direction e,-, i.e. 
the first component of the energy band i?( Q,k )'s Fourier trans- 
form, with the sum of quasi-momenta k extending over the 
first Brillouin zone of a lattice containing L sites. 

Note that the terms H e ,Ht and a part of 'Hj r do not couple 
different bands, whereas the on-site interaction term Hu\ oc 
conserves the local many-body parity 



II. TWO ALTERNATIVE DESCRIPTIONS 

We start with the single-particle Hamiltonian describing 
atoms of mass m in a 3D cubic optical lattice 
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with the same lattice depth s and spacing a in each dimension. 

We work in units of the recoil energy E r = Jr- (— J . Per- 
forming a band structure calculation and Fourier transforming 
the single-particle Bloch eigenstates leads to a multi-orbital ba- 
sis of Wannier orbitals [? ], for which we introduce the bosonic 
annihilation (creation) operators a, , Q (aj a ) at site i and in the 
band a = (a x ,a yi a z ). A short-ranged interaction for two 
atoms scattering in the s-wave channel only at the relevant 
energy scale can be well approximated by a <5-type contact 
interaction and characterized completely by the s-wave scat- 
tering length a s . The interaction strength parameter for the 
effective contact interaction is given by g = ATrha s /m and 
can also be expressed in the multiband Wannier basis as 
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where the matrix elements are defined in terms of the single- 
particle Wannier functions 

Uilt'X^ = §/ d 3 rw* uai (r)w* 2}a2 (r)w l3 , a3 (r)w l4 , ai (r) 

(3) 

Together with the contact interaction term, the full many-body 
interacting lattice Hamiltonian can be be written in terms of 
five contributions: 



(4) 



Ht + H-uXoc + n-t + 'Hu.nn + H\r 



Here, H e = Z) Q ( e(Q) _ A 4 ) J2i a l,a a i,a is tne on-site con- 
tribution of the single-particle lattice Hamiltonian Eq. {!} with 
e^ a ' being the mean energy of the band a and /i the chem- 
ical potential when switching to the grand canonical ensem- 

ble. The term H t = £ Q E( W )(<»U a + h - c is the 
tunneling between all pairs of nearest neighboring sites 
within the different bands a. = j.J2 k e iak e ' E (a M is 
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(for the x-dimension, others are analogous) along each dimen- 
sion, as shown in Appendix (A) . The local interacting Hamilto- 
nian at every site can thus be diagonalized in the subspace cor- 
responding to all multiorbital local states with the same parity 
as the ground state. The resulting eigenstates subsequently 
constitute an alternative set of basis states for the n-particle 
local Hilbert space. 

We now shortly recap the approximations made in the 
derivation of the standard Bose-Hubbard model: Firstly, one 
relies on a strong spatial localization of the single particle 
Wannier functions. For a sufficiently strong lattice depth s, 
this justifies to take only nearest-neighbor tunneling as well as 
only on-site interactions into account and neglect all others. 
Secondly, one assumes that all interband couplings (for any 
relevant operator) are negligible, thus justifying a truncation to 
the lowest single-particle band before constructing the many- 
particle Fock space. The presence of additional terms of the 
first kind is intrinsic to the problem and cannot be remedied. 
In a lattice of finite depth and with discrete translational sym- 
metry, the basis states in which all of these couplings would 
disappear necessarily has spatially completely delocalized ba- 
sis states (i.e. are the Bloch Fock states in the absence of inter- 
actions), which contradicts the initial goal of finding a spatially 
localized basis. On the other hand, the problem of interband 
couplings can be remedied by switching to the basis of local 
eigenstates. Here, the local interband couplings are contained 
to infinite order in each local eigenstate and the subsequent 
coupling between the different lattice sites gives rise to a new 
band structure, which we refer to as dressed bands. 

For a non-interacting system, the Fock state with \n)i par- 
ticles occupying the lowest Wannier orbital a = 1 at site i 
is the local lowest energy state and therefore well suited as 
a basis vector for a low-energy description of the system. In 
the presence of interactions |n)j is no longer the local lowest 
energy state, although the Fock states still provide a complete 
basis when taking all bands into account. However, for the 
simulation of interacting many-particle systems, one is often 
interested in a single-band description, and it is of great im- 
portance to find the best possible effective single-band basis. 
From these many-particle basis states one requires that they 

1. Have the highest possible spatial localization, i.e. mini- 
mize a localization measure such as the spatial variance 
of the density profile. 

2. Contain the local interaction-induced correlations, also 
lowering the many-particle energy expectation values of 
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FIG. 1. (Color online). Two-particle energy spectrum of the local 
Hamiltonian "n? in a 3D cubic 768nm lattice of depth s = 10 as 
a function of the interaction strength g. For any finite g the local 
ground state \ipo) ls no longer the pure Fock state, with all parti- 
cles occupying the lowest local Wannier orbital, but an admixture of 
higher Wannier orbitals is coupled by the local interaction terms. This 
lowers the total energy (including all orders in perturbation theory) 
and can be thought of as a dressed state in an effective lowest-band, 
as it evolves continuously from the g — limit and remains gapped 
from all higher dressed bands for typical interaction strengths. g va c 
is the two-body inter-atomic interaction strength of 87 Rb without the 
presence of an external magnetic field addressing the Feshbach reso- 
nance. In contrast to working in the truncated single orbital Wannier 
basis, this state gives a much better low-energy description contain- 
ing local correlations, as was confirmed experimentally 1 7 8 1. Spatial 
localization is furthermore guaranteed, since only higher Wannier or- 
bitals at the same site are occupied. 



these states (evaluated with the full interacting Hamil- 
tonian). 

3. Possess a well defined local particle number, such that 
the occupation number representation can be associated 
with these states. 

4. Are mutually orthogonal and span the complete low- 
energy subspace, i.e. formally constitute a basis. 

5. Recover the standard Bose-Hubbard model in the non- 
interacting limit. 

Having defined these requirements which we impose on an 
optimized effective single-band basis, the next task is to find a 
set of such states fulfilling the above requirements. We propose 
to use the many-particle eigenstates of the local interacting 
Hamiltonian T-Li oc — "H e + Hu^oc (which is a direct sum of lo- 
cal Hamiltonians H\ ), projected onto the Fock space spanned 
by the set of all non-interacting Wannier orbitals at a single 
site i. The above criteria are then fulfilled for the following 
reasons: 

1. Maximal spatial localization [? ] carries over from the 
maximum localization of the single-particle Wannier or- 
bitals at a given site. 

2. The multi-orbital, many-particle local ground state by 
definition minimizes the local energy and contains cor- 



relations in the interacting case, where the eigenstates 
are entangled with respect to the single-particle basis. 

3. Since the local truncated interacting Hamiltonian con- 
serves the local particle number 

KlT, a L% a ] = o, (6) 

OL 

(i) 

the eigenstates of H^' can all be chosen to have a fixed 
local particle number. For all typical experimental in- 
teraction strengths, the ground state is non-degenerate 
(thus necessarily possessing a fixed particle number). 
This allows a clear translation from the initial truncated 
single-band Wannier occupation into the new dressed 
band formalism: the initial local Fock state |rij)j is 
formally replaced by the local correlated ground state 
with n particles \ipo(n)) at the cost of renormalizing the 
Bose-Hubbard parameters. 

4. Local ground states with different particle number are 
orthogonal (or can be chosen as such in the case of de- 
generacy) since they are simultaneously eigenstates of 
the Hermitian local many-particle Hamiltonian. States 
at different sites on the other hand are orthogonal, 
since by construction they only occupy Wannier orbitals 
at different sites, which are orthogonal on the single- 
particle level. 

5. In the non-interacting limit, the local rt-particle ground 
state continuously converges to the local Wannier Fock 
state, thus recovering this limit. 

All longer-range matrix elements (beyond nearest neigh- 
bor) from both the lattice, as well as the interaction Hamil- 
tonian are contained in the long-range term Ti.\ r . These 
will not be discussed in further detail, since their transla- 
tion into the effective single band basis is identical to that 
of the nearest neighbor terms, they are however generally 
smaller in magnitude. The remaining terms (on-site and terms 
connecting nearest neighbors) from the interaction Hamilto- 
nian can be classified into four groups: on-site interaction 
terms forming the local interaction Hamiltonian %j/) oc = 

12 auaa , aa , at ,i U&^&s.tH. 4,a 1 a i, a2 a i,a 3 a i,a i ' the density- 
induced single-particle tunneling Hamiltonian H'{j nn contain- 
ing terms of the form a- „ a] „ a- „ a ■ „ + h.c. (as well as 
their counterparts under exchange i -H- j), pair tunneling 
terms of the form a\ ai al a2 a^ a3 a j ai + h.c. in 

nn' as 

well as the nearest-neighbor interaction Hamiltonian 

with terms a] „ a- „ at „ a„ „ + h.c. All terms of the latter 

*,CK2 J,c*3 J,c*4 

three types are contained in the nearest neighbor interaction 
Hamiltonian V lm , E Hu, rmV lawE = Hjj nn + Hij nn + H^ nn . 

III. DEFINITION OF THE LOW-ENERGY SUBSPACE 

In this section we systematically construct the effective low- 
energy subspace. The full Hamiltonian projected onto this sub- 
space gives the best possible description of interacting bosons 
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in a lattice at a sufficiently low temperature, where all higher 
dressed bands can be neglected. 

We diagonalize the local part of the Hamiltonian in Eq. 
which is a direct sum of local Hamiltonians 



fixed particle number and parity, leading to 

«£ = E^ ) ^m>«<V'S? ) | 



(8) 



Mum - E ^loc ( 



n 
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This can be achieved by diagonalizing the local Hamiltonian 
at each site separately, although for a homogeneous lattice the 
diagonalizations are of course identical and it suffices to per- 
form one only. Our formalism is however directly extensible to 
inhomogeneous systems, e.g. in the presence of additional spa- 
tial potentials or spatially dependent interactions. Note that in 

our notation operates only on the local Hilbert space of 
a single site (i.e. not on the complete lattice Fock space), the 
complete many-particle lattice Hilbert space is the direct prod- 
uct of all local Fock spaces over all sites. Truncating the local 
single particle space to the a max lowest bands, we diagonalize 
the local Hamiltonian in the Wannier Fock representation at 
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at site i: 
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FIG. 2. (Color online). Cartoon depicting our procedure to construct 
the localized low-energy, dressed band basis. The Wannier functions 
(shown for a lattice of depth s = 8E r ) are obtained from an ex- 
act single-particle band structure calculation on a large lattice, before 
one focuses on the diagonalization of the interacting many-particle 
Hamiltonian on the Fock subspace of multiorbital maximally local- 
ized Wannier states on a single site. This is not to be confused with a 
truncation of the spatial potential to a single lattice site. 



Due to number conservation of HyS , the eigenstates \tpm }i 
can always be chosen to be of fixed particle number n and 
contain an additional excitation index m. On a local level, 
a projection onto the many-particle low-energy space means 
only considering the local correlated local ground states 
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(9) 



ni,o,...,ni, 



with all different particle numbers n = n^o + • • • + "i,a m „- 
This low-energy projection can be extended to the whole sys- 
tem by defining the subspace spanned by the basis states 



m 



n i* 



(10) 



for all possible sets of integer local occupation numbers 
(m, . . . ,til). By construction, these states are all mutually 
orthogonal 



m ,n\ • • • w rtL ,n L i 



(11) 



which follows from the properties of the direct product in com- 
bination with the local states being different eigenstates of the 
same Hermitian Hamiltonian. It is also useful to define a low- 
energy projection operator 



v™= E i^" 1 ' 



(12) 



which projects any state from the full multi-orbital Fock space 
to the low-energy subspace of the full lattice. 



IV. TRANSFORMATION INTO THE NEW DRESSED BAND 
BASIS 

Having defined the effective single band space of interest, 
we now focus on expressing arbitrary operators in this sub- 
space. Here it proves very useful to define a set of new ladder 
operators 



k = E yftm 



(«)i 



(13) 



where i again refers to a physical site. It can be seen from 
the structure of these operators that any operator containing 
only transition elements between low-energy states of the type 
in Eq. |To) can be expressed in terms of these ladder operators 
and their Hermitian conjugates. Furthermore it can be directly 
verified that these operators fulfill bosonic commutation rela- 
tions 



b}]=S itj . 



(14) 
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Consequently, these ladder operators take over the role of the 
Wannier orbital creation and annihilation operators within 
a more appropriate single band description of an interacting 
bosonic lattice system. The next step is to express the original 
Hamiltonian and any other TV-particle operator in terms of the 
operators hi and b\, after projection onto the lowest dressed 
lowest band. This does not mean that all Wannier creation 
and annihilation operators Oj a and a\ are directly substi- 
tuted by hi and b\, but a systematic transformation is required. 
An arbitrary operator (acting on the full lattice) T> — Z)W 
expressed in terms of multiorbital lattice Wannier operators 
can be decomposed into normally ordered terms, where each 
term T>^> can contain operators corresponding to many dif- 
ferent lattice sites. Projecting this operator onto the lowest 
dressed band, i.e. multiplying with operator "P lowE from both 
the left and right, decouples this operator in the sense that the 
contribution to each lattice site can be considered individually. 
Omitting the site index, one such local term is thus generally 
of the normally ordered form 



A = 



(15) 



containing p creation and q annihilation operators and acting 
as the unit operator on all other sites. Thus, they transfer any 
basis state from the local (to + q — 1) -particle subspace to 
the (to + p — 1) -particle subspace for integer to > 1. We 
now use the fact that the set of all operators which change 
the local particle number of any state by some fixed integer 
number p — q and which annihilate all states containing less 
than q particles form a subspace of the complete operator 
vector space. We are however only interested in the low- 
energy subspace of this subspace, spanned by the operators 

{\4 p) )M q \ \4 p+1) u4 q+1) \, \4 p+ %M q+2 \ ■ ■ ■} 

4-p+m— 1 q-\-m— 1 

or equivalently by the set of operators {bl b i }. 

Thus we can express our local operator term as 



1owE U/ ^ ) q; 1 
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(16) 
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with coefficients h 
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, which are still to be de- 



termined. Here it proves useful to define the quantities 
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x(4 n+p -%l...aia....a \4 n) 



where we have dropped the site index i of which the f's are 
independent for a homogeneous lattice. The upper indices are 

labeled in a fashion, such that f; , , a a s is defined 

and can be non-zero for any integer r > 1. These coefficients 
can be directly calculated once the local eigenstates IV'o"'') are 
obtained from the exact diagonalization of the local Hamilto- 
nian. It is furthermore sufficient to restrict the indices to q > p, 
since all other cases are related by conjugation. For the spe- 
cial case of only annihilation operators being locally present 
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FIG. 3. (Color online). The coefficients fn as a function of inter- 
action strength g for a 738nm 3D cubic lattice. Note the symme- 
try relation f£l(2,o,o) = f«2(p,2,o) = /£=(o,o,2) and f^2(2,2,a) = 



/a=C2,0,2l = /a=(0,2,2) if the lattice iS isotropic. 



p(n) 



(i.e. p — 0), we define fjf^ ^ = /q p \ f° r notational 
convenience. These coefficients have the symmetry property 



(r) fr)* 
f(Pi,...,p,),( ai ,...,a p ) = f (ai,...,a p ),(fti,..., 



0, 



(18) 



and are furthermore invariant under permutations of indices 
within each bracket 

(r) (r) 

f(<Ti ai ),...,*(<*p)),(t(Pi),-MP<)) = J(oL 1 ,...,a P Up 1 ,...,p i y ( 19 ) 
Applying both sides of Eq. |l6) to the set of low-energy basis 
states | i/'o™'') for all integer n, using the definition in Eq. (17] 

and subsequently projecting onto (4 > o l+P |> we nn d the nn " 
ear transformation between the / and g coefficients 



(q 1 ,...,q p )(/3 1 ,...,/3,) 



V B h {m) 



with the matrix 



(n - 1)! 

Bn, m = 7 rr©(" - m) 

(n — my. 



(20) 



(21) 



and @{x) being the step function with 0(0) = 1. Explicitly, 
the first truncated part of the matrix is of the form 



B = 



( 1 
1 1 
1 2 2 
1 3 6 
1 4 12 






6 

24 








24 



1 5 20 60 120 120 



\ 



1 6 30 120 360 720 720 



V 



(22) 



The lower triangular form of the matrix is a consequence 
of the fact that the operator (6T) p+m (&) 9 + m annihilates any 
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state with less than (q + to) particles. In the non-interacting 
case, when the n-particle ground state is simply the state with 
all particles occupying the lowest Wannier orbital a = 0, 

we have /(al,...,a p )(ft,...A) = Sa ^° ■ ■ ■ Sa p,° S PuO ■ ■ ■ Sp t ,o 

, ap )(ft,...,ft) = S r,l S <* u ■ ■ ■ SapflSfofl . . . Sp qj0 , in 



local particle number operator transforms as 



(r) 



and ft ( ai ,... !ap )(ft,...A) - °r,l°ai,0 • ■ ■ °a p fiOp u Q ■ 

which case the density-induced transitions between the local 
ground states of different particle number vanish and the effec- 
tive low energy creation and annihilation operators are iden- 
tical to the Wannier creation and annihilation operators. 
To explicitly calculate the density-induced transition pa- 



rameters h, w fl 
which is found to be 
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we require the inverse matrix, 
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(24) 

It should be pointed out that due to this structure, the trun- 
cated inverse matrix is identical to the inverse truncated ma- 
trix. This finally allows us to express "Pi owE A'P lowE in terms of ef- 
fective low energy creation and annihilation operators, which 
explicitly reads 



^-/Sk + ^) (6t)P+2(&) 



9+2 



+ . . 



(25) 



with a. p = (ai, a.p) and f3 = (/3i, j3 q ). Note that in 
the non-interacting limit we have f R = 8 a ,o <5/3 ,o, all 

coefficients in brackets in Eq. [25] of higher order terms cancel 
I — 

and all operators ai >a i a [ a ) f° r tne l° wes t band can directly 
be replaced with bi (b\). 

For the specific case of a single annihilation operator 
'PiowEOa'PiowE the transformation is given by 



T^lowE a T^lowE 



(26) 



for which we show the first coefficients in Fig. as a 
function of the interaction strength. We point out that the 



r=l 



(6t)*-(6)*" = 6t 6j (27) 



since ^ Q f} a {r a \ = 1 for all r and the property of the trans- 
formation matrix in Eq. j^j ~^2 n ( B 1 ) m ,n — $m,i- Hence the 
local particle operator in our single dressed band counts all 
particle in all local orbitals. 

The transformation relations are the central result of this 
section: given any physical operator in the single-particle 
Wannier basis (as is usually the case) which acts on the sys- 
tem or is measured, one cannot simply substitute the Wan- 
nier creation and annihilation operators ai M with the bi op- 
erators from the effective single-band model with density de- 
pendent parameters! Rather, the transformation relations from 
Eqns. jl6|25} have to be used to systematically express this op- 
erator in the effective low-energy subspace. 

In the following sections we will apply this transformation 
to the various terms of the original Hamiltonian. The purely 
local terms exactly recover the effective multibody interac- 
tions found in II 2 3 II and are diagonal in the dressed band ba- 
sis with density-dependent interaction parameters. Non-local 
nearest neighbor terms originate from Tit + Hu,rm and are 
non-diagonal in the new basis, leading to the low-energy rep- 
resentation 



' lowE Ti tot 'lowE 'lowE Ti\oc ^lowE nn 



njint 
n U,rm 



(28) 



of the full initial Hamiltonian. We neglect all long-range (be- 
yond nearest neighbor) terms contained in Ti\ x in this work, 
their transformation is however identical to that of the near- 
est neighbor terms. Thereafter we will discuss an equivalent 
formulation of the extension of the Bose-Hubbard model using 
density-dependent parameters for the various terms at the cost 
of additionally summing over the set of all local low-energy 
states. 



V. APPLICATION TO THE BOSE-HUBBARD MODEL: 
MULTI-BODY INDUCED TUNNELING AND INTERACTIONS 

In this section we discuss the four relevant terms in the 
dressed single-band Bose-Hubbard model. This contains all 
relevant local and nearest neighbor processes. The amplitudes 
for processes on neighboring sites in the multi-body induced 
picture are shown in Fig. 



A. Single-particle tunneling term 

We now have to gather all operator contributions in the 
total Hamiltonian Ti to , that give rise to single-particle tun- 
neling transitions between nearest neighboring sites i and j. 
Clearly, the tunneling term from the original single-particle 
lattice Hamiltonian Tit is such a term, but also the interaction 
Hamiltonian Ti.jj.xm contains single-particle transition terms of 
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FIG. 4. (Color online). Amplitudes for all multibody-induced transitions on nearest neighboring sites. The lowest order, most relevant processes 
are shown as functions of the lattice depth s at a fixed interaction strength g = 5g vac = 0.186i5 r a 3 for 87 Rb in the upper row, as well as their 
dependence on the interaction strength g at a fixed lattice depth s = 10E r in the lower row. Note the different scaling of individual graphs 
for visual clarity, which is given in the legends. The plots (al,a2) in the left column show the effective single-particle tunneling, the dotted 
and solid lines show the contributions from the interaction term Hu, nn and the single particle particle tunneling term T-Lt respectively. With 
increasing order (mi, m,2) the contribution from the interaction becomes more important and eventually dominates. Whereas the dependence 
on the interaction strength g is monotonic, the dependence on the lattice depth is considerably more complicated, being non-monotonic and 
even leading to sign changes. We point out that as a function of the lattice depth s, the magnitude of the higher order terms is most significant 
in the region of s ~ 10E r or slightly below, which is also the relevant region for the superfluid-Mott insulator transition (depending on the 
interaction strength g, i.e. the Feshbach resonance). In the noninteracting limit (7 — ^ all terms, except the lowest order single particle tunneling 
Mi 1 vanish, recovering the usual lowest band Bose-Hubbard model. However, at any finite interaction strength with other terms becoming 
non-zero, the Bose-Hubbard model truncated to the lowest single particle Bloch band does not give the correct low-energy description. 



this type, which we denote by Hjj nn . For each fixed set of 
nearest neighbor sites i and j there are four relevant terms for 
this process and thus the total single-particle hopping Hamil- 
tonian originating from the interaction term is 



*3, CK4 w i,ftri,ft2 Lt !,ft3 li j,ft4 



iTj(i,jJ,j) J J a „ 1 + V, C 

(29) 

Note that in contrast to Ti. t , which is diagonal in the band in- 
dex, T-Lfj nn couples Wannier states in different bands on neigh- 
boring sites. 

The single-particle tunneling %j can also directly be writ- 
ten in terms of 6-operators, giving rise to multiparticle induced 
single-particle tunneling 



with the (s, r)-particle-induced tunneling amplitude 



M, 



= E^ S) ^" ) + E 



Tj(i,i,i,j) 

Oil ,C*2 ,CK3 ,CK4 



x h^^ 1 ^" |,W 1 ff('JJj) 



Ql ,C*2 ,C*3 ,C*4 



(31) 



(a 3 )(a 2 ai) on 



ai ,a 2 ,03,04 (a 2 )(a 3 («4) ai 



For r = s = 1 Eq. poj is simply a usual single-particle tun- 
neling term, containing all multi-orbital contributions in "Ht, 
but no contribution from %'},„„, since h, , vanishes 

for any r < 1. However, there are also additional multibody- 
induced single-particle tunneling terms present: for r > 1 or 
s > 1 a single particle can tunnel between neighboring lattice 
sites with an amplitude M r s if r — 1 and s — 1 additional 
particles (additional to the one tunneling) are present on the 
lattice sites. We therefore refer to these processes as being 
multibody-induced. 



00 

r,s=l 



(30) 



B. Two-particle correlated hopping 

In contrast to the non-interacting lattice Hamiltonian, the 
interaction term "H[/,nn also contains two-particle correlated 
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tunneling transition elements in the term H(j nn - A single ap- with e a the mean energy of the ath Bloch band and 



plication of such a term to the state l^o™*^)*!^ )j leads to a 
correlated tunneling of two particles on neighboring sites i and 

j, leading to states of the form |^q" * +2 '*)i|'0o™ J Clearly 
such operator terms are beyond the standard Bose-Hubbard 
model and cannot be contained in a renormalized tunneling 
parameter. They may however lead to interesting effects and 
we additionally include them in an extended description of the 
interacting lattice model. 



£r — E £ah, a 



(r) 

(COW 



(37) 



Note that 



a a 



(38) 



eo, 



hs+ifh \s-i/7.1v-i^ .) r + x ] +h.c. 



(32) 



with the (s, r)-particle-induced two-particle tunneling ampli- 
tude 



H — TjihhJj) /,<>)* hi 

^s,r — 2_^, u a 1 ,a 2 ,a 3 ,a 4 ,"'a 1 ,a 2 "'c 



"3, "4- (33) 



Ql,Q2 ,C*3 ;C*4 



C. Nearest neighbor interactions 

Counting terms corresponding to nearest neighbor inter- 
actions in the total nearest neighboring part of the Hamilto- 

„ „ „ „ UX-, ao a, ua Oj ~ €Li ™ ™ ™ With 

CKl ,Q2 1<^3 -Ct4 1 ' 2 ' 3 ' 4 *2i<^2 *3,CK3 ^4 5 C*4 

all i m being of two nearest neighbor sites i and j, there 
are four terms corresponding to nearest neighbor interactions. 
These are all equivalent and after permuting indices the full 
nearest neighbor interaction Hamiltonian can be written in the 
form 

njint _ a \ " \ " Tj{i,i,3,j) 
IL U,rm / , / j ^ a 1 ,a 2 ,a 3 ,a 4 



«l,«2,a3,(«4 



(34) 



r,«=l {id) 



where the last line is in the effective multibody nearest neigh- 
bor interaction picture with the parameters 



since the ground state containing one particle \ipQ ) is not af- 
fected by interactions and thus is the single particle Wannier 
state. The on-site interaction term expressed in terms of b op- 
erators is of the form 

oo 

V^H UM V^ = E E U r (bl) r+1 (bi) r+1 (39) 



with 



U = V r/frw) h {r) 

Ur ~ Z_^i U ai,a 2 ,a 3 ,a i n (a 1 ,a 



2)(a3.«4 



(40) 



Oil ,Cf2 ,Q!3,Q;4 



Both on-site terms can be combined to obtain the local part of 
the Hamiltonian in terms of the b's 

oo 

V^n uM V lowE = E(eo - fi)b\lH + E E V r(4) r (h) r 



r=2 



(41) 

with V r = e r + U r -\. The V's are related to the eigenenergies 
of the low-energy basis states \ipQ^), sm ce 

v^ s u loc v lms -EE 4 n) \4 n) )u(4 n) \- («) 



Letting both sides of Eq. p?| act on the local low energy basis 

states \ipQ )i f° r a U integer n leads to the expression of the 
higher order interaction amplitudes in terms of on-site many- 
particle eigenenergies 



E, 



v rn — / , \*- > )m,n 

' 71 



(n) 



(43) 



VV Str -1 U *ua 2 , a3 , ai n ( ai )( a2 ) n ( a3 )( a4 )- (35) 



Ql ,Q2 .,Ct3 ,C«4 



VI. DENSITY-DEPENDENT PARAMETER FORMULATION OF 
THE BOSE-HUBBARD MODEL 



D. On-site terms 
The single-particle on-site terms are given by 

(36) 

= £(eo-M)&&+£X>$)W" 

i i r—2 



An alternative and equivalent description to the multibody 
induced picture is a formulation in terms of density-dependent 
parameters. In this picture, one performs a summation over 
all local occupation numbers and allows the matrix elements, 

which are coefficients of operators of the form |V'o i )«(' ! /'o i I 
to acquire a density dependence beyond the bosonic statistical 
factor. This picture is convenient to directly infer certain mul- 
tiorbital effects, such as the observed energy peak positions 
in quantum phase revival spectroscopy experiments [7] HI- 
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However, the density-dependent representation is not always 
the most convenient approach for treating multiorbital effects 
with usual many-body methods. For instance, bosonization 
or even the site decoupling mean-field theory cannot be per- 
formed within this framework. The effective representation 
discussed in the previous section Sec. (v) , which does not re- 
quire an external summation over all local occupation num- 
bers is more favorable in this sense and has been successfully 
applied for the pure on-site terms in describing effective multi- 
body interactions II23II and commonly accepted as the repre- 
sentation of choice [7] [H . 



A. Single-particle tunneling term 

We first seek the single-particle tunneling term in the 
Hamiltonian with a density-dependent tunneling parameter, 
i.e. of the form 

(i,j)ni,nj (44) 

®^n j \^- 1) ) jj {^ ) \+h.c. 

In this section we omit writing the unit operator on other sites 
for any local operator: i.e. any operator A^> acting only on the 
local Fock space of site i is to be implicitly understood as being 
extended to the full lattice Fock space as A^' <8> n®j=£i ^j- 
Each operator term in Eq. {44} can also be written with the 

use of local projectors Vn = |V , o"^)«(' , /'o™'' I on the n-particle 



ground state at site i as 

The final summation over all occupation numbers rij , rij is 
however always necessary in the density-dependent param- 
eter representation of any operator. 

To transform the total single-particle tunneling Hamiltonian 
%j = H t +Hy nn to the density-dependent parameter form of 
Eq. (44) , we insert the unit operator of the low-energy subspace 
on both the left and the right. For one given operator term, 
such as Ai.j = a\ ai al a2 a i a3 a j ai , the expectation reduces 
to a product of expectation values at different sites 

(46) 

where the coefficients f$ n, \ r \ are defined in Eq. Jl7b. 
On all other lattice sites different from i or j, Aij acts as 
the unit operator. The same procedure can be used for %f 
We furthermore use the fact that for a time reversal sym- 
metric system, all Wannier functions can be chosen purely 

real and the matrix element Ua^'a 2 2 faC,a 4 is invariant under 
the 24 possible permutations of index pairs (i n ,a n ). To- 
gether with the discrete translational symmetry we thus have 
Uui,aiia 3 ,a4, = Ua^a2,a3,a! and upon relabeling the sum- 
mation indices ai W 04 obtain the total density dependent 
single-particle tunneling parameter after collecting all terms 
ofHj 



I - P 



jU = .(a) f ( nj ) f(n 



E 

a 1 ,Q2,a 3 ,a 4 



Ot\ ,a 2 ,Ct3 ,Q4 



A n i)* 

'■ J (a 3 )(a2,ai) J 



'■> (ot 2 )(a 3 , 



(47) 



B. Two-particle correlated hopping cients 



For a density-dependent representation of the two-particle ni,rlj _ „ ai,a2 ' a3 '° 4 ai '" 2 Jq3 '" 4 ( 49 ) 

nneling parameter, they art 
term in the total Hamiltonian 



tunneling parameter, they are represented by an additional 



C. Nearest neighbor interactions 

^a,nn = E E ^,%VK + l)(ni + 2) |^ n * +2) )«(4 n,) i 

n t ,n 3 jj ie neares t neighbor interaction Hamiltonian in its density- 

® ^Ujinj-l) |Vo"^ 2) )« I + h - c - dependent parameter representation reads 

(48) -EE ^ * \^ ] )M nih 



(i,i)«i,%- (50) 

The same procedure as described in Sec. jV A\ can be used r n \ r n .\ 

to obtain the density-dependent two-particle tunneling coefR- 
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Note that the coefficient / 



(n) 

(01,02) 



in this case is exactly the 



local single-particle density matrix in the multi-orbital Wan- 
nier representation for the local many-particle ground state 

\4 n) )- 



VII. ENERGY REDUCTION MECHANISMS: COMPARISON 
BETWEEN ORBITAL DEFORMATION AND HIGHER ORDER 
CORRELATION EFFECTS 

We demonstrate now that the deformation of single-particle 
orbitals by the interactions is not the main effect to lower the 
local on-site energy. The higher order correlations, which can- 
not be understood as such a deformation, are the dominant ef- 
fect to lower the energy. Therefore, a single-particle picture 
and wave functions are not sufficient for understanding the 
effect of interactions on the local level, since entanglement be- 
comes important. This can best be seen in the two-particle 
correlation function in Fig. {6}. 

Given the state \ip^)i, the local single-particle density ma- 
trix is hermitian and can thus be expressed in terms of its or- 
thogonal eigenvectors (corresponding to single-particle states 



in the respective basis) 
tive eigenvalues A; 



and the corresponding real, posi- 



(n)\ 



(51) 



The eigenvalues A; and the associated single-particle states do 
of course not depend on the basis in which the single-particle 
density matrix is evaluated in. 

We can now construct an artificial state 



(52) 



for comparison, which leads to the identical single-particle 
density matrix, but does not contain the higher order corre- 
lations. 

Here we defined the creation operators for the eigenstates of 
the single-particle density matrix dj = ^ Q 0« aj . The state 
in Eq. |52) can be thought of as having the identical single- 
particle properties as the true local ground state, and would 
be the most natural many-particle state for thinking in terms 
of spatially broadened single-particle orbitals due to the in- 
teractions, as commonly referred to in literature [7j. It does 
however not contain the same higher order correlations as the 
original state, for instance the two-particle correlation func- 
tion G^, = (ip uc \d] d\,d l ,d l \ip uc ) = Xi Sim for n — 2 parti- 
cles or more. In Fig. {5} the energy expectation value of the un- 
corrected state E uc = (^puc\H\' l Puc) is compared to the true 
ground state energy as a function of the interaction strength 
g. Since the above constructed state \ip U c) becomes the true 
ground state in the non-interacting limit, the energies agree 
in this limit. However, significant deviations arise at finite 
interaction strengths g, indicating that the simple picture of 




0.1 

(in units of E r a ) 



FIG. 5. (Color online) Comparison of the true many-body local 
ground state energy E gs and the energy E uc of the artificially cre- 
ated state \ ipuc) of Eq. (52) with the same single-particle density ma- 
trix (i.e. the same broadened orbitals and their occupation), but no 
higher order correlations. The deviation between the two becomes 
very significant with increasing interaction strength and is well in the 
experimentally observable regime (here for s = 10, but the effect be- 
comes stronger with increasing lattice depth), showing that the usual 
simple picture of broadened single-particle orbitals is insufficient to 
explain the effects of interactions on the local many-body state. 



spatially broadened single-particle orbitals cannot explain the 
main energy reduction mechanism. 

The significant change of the on-site many-body state lies 
in the higher order correlation functions, with particles mutu- 
ally reducing their spatial overlap. It is not contained in and 
cannot be understood on single-particle level (since all single- 
fa) 

particle properties of the local state are contained in p a ' a )- To 
substantiate this point, we calculated the normalized second- 
order correlation function 



<? (2) (r,r') = 



(53) 



(2) 

for a local interacting two-particle ground state |"0q ), shown 
in Fig. |6j. The normalized second order correlation function 
can be understood as a conditional probability: for the non- 
interacting state \n = 2) we have </ 2 '(r, r') = 1 — ^ = \, 
which is our reference and which we refer to as uncorrelated 
by interactions. A value of g^ 2 ' (r, r') < 1 indicates a reduced 
probability for a particle to be found at location r' if another 
particle is located at r or vice versa and is therefore anticorre- 
lated in this sense. This anticorrelation can be seen along the 
diagonal line x — X 1 in Fig. |6j, where two repulsively interact- 
ing atoms have a reduced probability to be found at the same 
spatial position x = x' . Since all particles are restricted to oc- 
cupy Wannier orbitals at the same site, the conditional proba- 
bility g( 2 )(r, r') has to be increased elsewhere, i.e. correlated. 
With increasing interactions, the specific shape of </ 2 )(r,r') 
is independent of the interaction strength g in this regime in 
the sense that the deviation from the uncorrelated case scales 
linearly with g. This can clearly be seen by comparing the 
different functions in the left column of Fig. {6}. In contrast, 
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g (2) (x,x') G (1) (x,x') 




FIG. 6. (Color online). Correlation functions g^ (x,x') (left column) 
and single-particle density matrix G' 1 ' (x, x') (right column) in real 
space at y — y' — z — z' = for a single site in units of the lattice 
constant a. Results are shown for n = 2 87 Rb atoms in a 768nm 
s = 10£V lattice for three different interaction strengths, up to four 
times the vacuum interaction strength g m . The lower graphs are cuts 
along the dotted lines in the plots above. With increasing interaction 
strength, the atoms mutually avoid each other, as can be seen in the 
decrease of g 2 \x,x) along the diagonal (x = x') and an increase 
for [a; — x | ^ 0.2, as compared to the non-interacting case where 
g^ 2 \x, x') = 0.5 for n = 2. The effect on the density distribution 
along one direction is significantly weaker, with the main effect being 
a reduction of density at the center of the lattice site. 



VIII. CONCLUSIONS 

We have formulated a systematic derivation of an effective 
low-energy, single-band basis for ultracold bosonic atoms in 
optical lattices in the presence of interactions. Some prop- 
erties intrinsic to our formalism, such as density-dependent 
interaction parameters or the appearance of effective multi- 
body interactions have been previously discussed and experi- 
mentally confirmed. We introduce ladder operators fulfilling 
bosonic commutation relations within the new dressed band 
basis, which are shown to be the bosonic operators used within 
an effective Bose-Hubbard model for the system. It is however 
shown that these are not the original lowest band Wannier cre- 
ation and annihilation operators beyond lowest order and we 
derive a simple prescription for the transformation of arbitrary 
operators into the new low-energy dressed band basis. These 
transformations are used to systematically treat all terms in 
the interacting lattice Hamiltonian and give rise to multibody- 
induced single and pair particle tunneling, as well as multi- 
body local and nearest neighbor interactions. The amplitudes 
for these processes are calculated and compared to renormal- 
ized parameters in the density-dependent representation of the 
Bose-Hubbard model. The latter formulation, although fully 
equivalent, is however less favorable for the treatment with a 
number of common theoretical methods, since it contains an 
external summation over the set of all n-particle states at each 
lattice site. We furthermore show that the commonly used 
single-particle picture of spatially broadened Wannier orbitals 
cannot describe the observed energy reduction of the local 
many-body state. The relevant mechanism is mutual avoid- 
ance of the various atoms at a given lattice site, which is a 
many-particle effect contained only in the higher order corre- 
lation functions. 

Note added: Related results on the calculation of density- 
dependent single-particle tunneling parameters have recently 
been reported in Ref. [24] during the final stage of this work. 
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Appendix A: Parity Symmetry 



the density profile shown in the right column is only slightly 
changed by the interactions. The main effect is a reduction of 
the maximal density at the center of the lattice site, whereas 
only a minimal broadening of the density profile is visible. 



Here we prove that the local interacting Hamiltonian pre- 
serves the multiparticle parity along each direction at a given 
site i and drop the site index for this section. The local Wan- 
nier orbitals are either fully spatially symmetric or antisym- 
metric in all of the three spatial dimensions, thus render- 
ing the multiparticle parity operator diagonal in the Wannier 
Fock space representation. Since the local lattice Hamiltonian 
is also diagonal in this representation, we can directly infer 

VU,qW] =o. 
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We now focus on proving the second relation 
[^(7,loc) Q^] = 0. This is equivalent to the statement 
that both operators share a common basis of eigenvectors, or 
equivalently, that in a given basis one operator is block diag- 
onal, with non-diagonal blocks only within those subspaces 
where the other operator is a scalar multiple of the identity. 
We take the basis of Wannier Fock states and consider the 
expression 

(H UM Q {X) - Q {x) H uM )\n ai ,. . .,n aM ). (Al) 

Since the state is an eigenstate of Q^ x \ the first term corre- 
sponds to XHu\oc\ n ot\, ■ ■ ■ i n a M ), where 



£« a ,=i,3,...£c 



A = (-l) 



is the eigenvalue belonging specifically to this state. 



(A2) 



We now consider the second term 
2 



' n U ai ,a 2 ,a 3 ,a i a>a i a>a 2 a- a3 a a4 \na 1 , ■ ■ ■ T U a M ) 

Ql ,CK2 ,CK3 ,Q4 

(A3) 

and use the property that the interaction matrix element 

U ai ,a2,ot3,a!4 — T\d=x,y,zUa 1 , d ,a 2 ,d,a3,d,a4,d factorizes into a 
product of terms from the individual dimensions. One such 
term 



(A4) 

vanishes if it contains an odd number of odd functions. Conse- 
quently all non-vanishing states created in the sum in Eq. jA3| 
are also eigenstates of to the same eigenvalue A as the 
initial state. 

As th e Bloch Fock states constitute a complete basis, 
Eq. i A3 i holds on an operator level, i.e. [Hi/\ oc , Qv*'] = 0. 



[1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 
(2008). 

[2] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, 

Phys. Rev. Lett. 81, 3108 (1998). 
[3] R. Feynman, Int. J. Theor. Phys. 21, 467 (1982). 
[4] W. Hofstetter, J. I. Cirac, P. Zoller, E. Dernier, and M. D. Lukin, 

Phys. Rev. Lett. 89, 220407 (2002). 
[5] T. Esslinger, Annu. Rev. Condens. Matter Phys. 1, 129 (2010). 
[6] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansen, and I. Bloch, 

Nature 415, 39 (2002). 
[7] S. Will, T. Best, U. Schneider, L. Hackermiiller, D.-S. Luhmann, 

and I. Bloch, Nature 465, 197 (2010). 
[8] S. Will, T. Best, S. Braun, U. Schneider, and I. Bloch, Phys. Rev. 

Lett. 106, 115305 (2011). 
[9] S. Trotzky, L. Pollet, F. Gerbier, U. Schnorrberger, I. Bloch, N. 

V. Prokofev, B. Svistunov, and M. Troyer, Nature Physics 6, 998 

(2010). 

[10] S. Ospelkaus, C. Ospelkaus, O. Wille, M. Succo, P. Ernst, K. Sen- 
gstock, and K. Bongs, Phys. Rev. Lett. 96, 180403 (2006). 

[11] T. Best, S. Will, U. Schneider, L. Hackermiiller, D. van Oosten, I. 
Bloch, and D.-S. Luhmann, Phys. Rev. Lett. 102, 030408 (2009). 



[12] M. J. Mark, E. Haller, K. Lauber, J. G. Danzl, A. J. Daley, and 

H.-C. Nagerl, |arXiv:l 107. 1803^ 1. 
[13] W. S. Bakr, P. M. Preiss, M. E. Tai, R. Ma, J. Simon, M. Greiner, 

larXiv:1105.5834H . 
[14] J. Heinze, S. Gotze, J. S. Krauser, B. Hundt, N. Flaschner, D.-S. 

Luhmann, C. Becker, K. Sengstock, |arXiv: 1107.2322| fl. 
[15] D.-S. Luhmann, K. Bongs, K. Sengstock and D. Pfannkuche, 

Phys. Rev. Lett. 101, 050402 (2008). 
[16] R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev. A 79, 

011606 (2009). 

[17] A. Mering and M. Fleischhauer, Phys. Rev. A 83, 063630 (2011). 
[18] M. Cramer, S. Ospelkaus, C. Ospelkaus, K. Bongs, K. Sengstock, 

and J. Eisert, Phys. Rev. Lett. 100, 140409 (2008). 
[19] M. Cramer, Phys. Rev. Lett. 106, 215302 (2011). 
[20] M. Snoek, I. Titvinidze, I. Bloch, and W. Hofstetter, Phys. Rev. 

Lett. 106, 155301 (2011). 
[21] K. R. A. Hazzard and E. J. Mueller, Phys. Rev. A 81, 031602(R) 

(2010). 

[22] O. Dutta, A. Eckardt, P. Hauke, B. Malomed, and M. Lewenstein, 

New J. Phys. 13, 023019 (2011). 
[23] P. R. Johnson, E. Tiesinga, J. V. Porto, and C. J. Williams, New J. 

Phys. 11, 093022 (2009). 
[24] D.-S. Luhmann, O. Jiirgensen, and K. Sengstock, 

larXiv:iT0 8.3013yl. 



